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ABSTRACT 


The  paper  is  concerned  with  a  theoretical  study  of  shearing  flow 
bounded  by  a  wavy  surface  with  consideration  of  turbulent  flow  above 
the  air-sea  interface.    Account  is  to  be  taken  of  turbulence  in  the  flow 
through  the  use  of  Reynolds  stresses  associated  with  turbulent  flow. 
An  adaptation  of  the  mixing  length  theory  as  applied  to  pipe  flow  is  made 
for  channel  flow  and  the  resulting  mixing  length  versus  energy  relation- 
ship is  incorporated  in  the  Reynolds  stress  term.     Finally,  the  rate  of 
wave  growth  is  calculated  from  the  normal  surface  pressure  in  phase 
with  the  wave  slope. 

A  curvilinear  coordinate  system  which  follows  with  the  wave  train 
is  used  in  order  to  simplify  the  formulation  of  the  problem.    All  para- 
meters are  non-dimensionalized  and  the  analysis  is  made  considering 
a  velocity  profile  adapted  from  pipe  to  channel  flow. 
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SYMBOLS  AND  ABBREVIATIONS 

a  non-dimensional  wave  amplitude 

A  constant  parameter 

c  non-dimensional  wave  celerity 

C  constant 

d  position  at  which  upper  boundary  conditions  are  applied 

E'  turbulent  energy 

E  wave  energy 

E  rate  of  energy  input  to  wave 

F  function  determined  by  solution 

Froude  number 


r 


g  acceleration  of  gravity  (980  cm  sec     ) 

H  channel  height 

J  Jacobian 

K  kinematic  eddy  viscosity 

mixing  length 

mixing  length  (perturbation  term) 

mixing  length  (zero-order  term) 
P  pressure 

q  velocity  of  main  flow  above  boundary  layer 

U  velocity  of  main  flow  in  the  channel 

u  non-dimensional  velocity  parallel  to  \     axis 

u*  friction  velocity 


SI 

2, 

L 


v  non-dimensional  velocity  parallel  to  V\    axis 

°<^  constant 

X  real  component  of  pressure 

o  imaginary  component  of  pressure 

€  real  component  of  shear  stress 

£  imaginary  component  of  shear  stress 

o  constant  parameter 

X  mixing  length  constant  (0.4) 

A  non-dimensional  wave  length 

^  stream  function 

AU.  molecular  viscosity 

B^  air  density  (0.0013  g  cm"^) 

JJj  water  density  (1.027  g  cm-3) 

L  shear  stress 

1*  7)  curvilinear  coordinates 
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1.     INTRODUCTION 

A  semi-empirical  study  of  turbulent  air  flow  over  a  wavy  surface 
is  made  in  which  normal  and  tangential  stresses  on  the  surface  are 
examined  with  the  object  of  estimating  the  growth  of  water  waves  due 
to  a  turbulent  wind . 

Ursell  (195  6)  gives  a  summary  and  a  critical  review  of  the  theories 
of  wind-wave  formation  extant  in  1956.    The  review  gave  the  necessary 
impetus  for  further  investigations,  both  theoretical  and  experimental, 
in  this  area.    Since  1956  there  has  been  much  work  done  on  the  generation 
of  water  waves  by  surface  stresses  induced  by  the  wind.     Miles  (195  7) 
and  Benjamin  (195  9)  considered  laminar  air  flow  over  a  wavy  water 
surface  and  deduced  formulae  for  exponential  wave  growth  due  to  an 
instability  mechanism.     Later  Phillips  (1966)  included  air-stream 
turbulence  in  the  instability  model  from  a  semi-empirical  point  of  view. 

In  his  initial  work,  Miles  (1957)  assumes  an  inviscid,  incompress- 
ible, two-dimensional,  parallel  shear  flow  coupled  with  a  prescribed 
two-dimensional  deep-water  gravity  wave.     Turbulent  fluctuations  in 
the  perturbation  equations  of  motion  are  neglected,  even  though  they 
are  decisive  in  maintaining  the  mean  shear  flow.     Furthermore,  the  wave- 
induced  velocities  and  pressure  are  small,  which  justifies  linearization 
of  the  equations.    The  water  is  also  assumed  to  be  inviscid,  incompress- 
ible, and  irrotational,  and  the  slope  of  the  displaced  surface  is  assumed 
to  be  small.     Only  that  component  of  the  aerodynamic  force  in  phase 
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with  the  wave  slope  is  important  in  the  wave-generation  mechanism. 
The  viscous  drag  forces  of  the  air  are  assumed  to  be  negligible  compared 
with  the  normal  pressure  in  their  effect  on  the  surface  wave. 

This  model  is  used  to  determine  a  first-order  approximation  to  the 
disturbed  motion  of  the  air  and  the  consequent  energy  transfer  to  the 
wave.    The  velocity  distribution  is  taken  to  be  that  for  turbulent  flow 
over  an  aerodynamically  smooth  water  surface.    The  result  rests  on  the 
solution  of  the  inviscid  Orr- Sommerfeld  equation  which  has  a  singularity 
at  the  height  where    \J       =    c      (the  critical  level) .    An  exponential 
increase  in  amplitude  of  water  waves  with  phase  velocity    c    occurs  if 
the  mean  air  flow  velocity     U  (-y)        in  the  direction  of  wave  travel  varies 
with  height  above  the  surface  such  that  U       at  the  critical  level  is 
negative.     This  is  almost  always  the  case. 

Phillips  (1966)  attempts  to  include  wave-induced  turbulence  in  the 
instability  model  by  including  in  the  momentum  flux  equation  an 
appropriate  integral  over  the  entire  layer  of  air  of  the  contribution  due 
to  perturbation  eddy  viscosities.    This  integral  is  multiplied  by  an 
unknown  constant    A    which  is  estimated  from  experimental  results  from 
flow  over  fixed  wavy  surfaces.    The  uncertainty  in    A    is  estimated  to 
be  +  50  per  cent.     Phillips  concludes  that  even  though  there  exists 
a  high  degree  of  uncertainty  in    A    the  results  do  provide  a  framework 
in  which  experimental  measurements  may  be  interpreted. 

According  to  Miles  (1967)  ,  experimental  results  indicate  that  a 
theoretical  model  based  on  laminar  flow  may  be  adequate  in  the 
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laboratory  but  not  on  an  oceanographic  scale.    The  inviscid,  laminar  model 
underestimates  the  energy  transfer  from  wind  to  waves  over  at  least  a 
significant  portion  of  the  spectrum  for  an  open  sea.    This  suggests  that 
the  significance  of  wave-induced  perturbations  in  the  turbulent  Reynolds 
stresses  for  momentum  and  energy  transfer  from  wind  to  wave  must 
increase  with  an  appropriate  scale  factor.     Miles  concludes  that  further 
theoretical  work  in  this  area  appears  to  demand  some  'ad  hoc'  hypothesis 
to  describe  these  perturbation  Reynolds  stresses. 

This  can  be  done  by  considering  a  slightly  generalized  form  of  the 
classical  mixing-length  hypothesis.    A  means  must  be  found  to  predict 
the  value  for  the  change  in  mixing  length  upon  which  the  perturbation 
Reynolds  stresses  are  based.    Zagustin  (1963)  has  given  a  semi-empirical 
solution  for  turbulent  flow  in  a  pipe  which  gives  good  agreement  with 
experimental  data  and  which  predicts  the  value  of  the  mixing  length. 
This  model  is  generalized  below  to  include  flow  over  a  wavy  surface  in 
a  manner  analogous  to  that  of  Miles  and  Benjamin  in  that  it  assumes 
the  flow  to  follow  the  wavy  surface.    Viscous  stress  terms  are  replaced 
by  the  Reynolds  stresses  associated  with  turbulent  flow.    Normal  Reynolds 
stresses  are  neglected.    Zagustin' s  method  of  closing  the  system  of  mean 
turbulent  equations  of  motion  is  not  unique  (see,  for  example,  Squire, 
195  9) ,  but  there  seems  to  be  no  good  reason  to  use  any  other  model  in 
its  place.    The  investigation  below  details  the  application  of  Zagustin' s 
work  to  this  model  and  outlines  the  solution  of  the  differential  equations 
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using  numerical  methods  and  appropriate  boundary  conditions.    The 
unknown  parameters  in  the  model  are  fitted  to  experimental  data  describing 
the  normal  pressure  in  phase  with  the  wave  slope. 
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2.     FORMULATION  OF  THE  MODEL 

The  zero-order  flow  is  taken  as  flow  in  a  two-dimensional  channel. 
All  lengths  are  made  non-dimensional  using  the  channel  halfwidth,     H. 
Velocities  are  made  non-dimensional  using  the  midstream  velocity,  Uc 


All  stresses  are  made  non-dimensional  by  dividing  by    o  \J0      ,  where  o 
is  the  density  of  air. 


2.1    THE  ZERO-ORDER  FLOW 

The  velocity  is  expressed  in  terms  of  a  stationary  and  a  fluctuating 
(turbulent)  part: 


(2.1)  ttJL       =     "0"       ■+■      WL 


Following  Zagustin,  the  Reynolds  shear  stresses    L     are  written  as 

2. 


(2.2) 


"f     =     L     (qrcLcLq) 


where    Q      is  the  mean  speed.    This  can  be  considered  a  defining 
equation  for    L       (although  here  it  appears  as  the  usual  mixing  length) 
in  that    L     will  not  be  prescribed,   but  will  be  part  of  the  solution.    Then 
the  two-dimensional  mean  equation  of  motion  is 

(2.3)  ( "0"    •  grad)  \J         =    -grad  P    +  |K  x  grad  T 

A 

where     IK       is  a  unit  vector  normal  to  the  plane  of  the  flow,  P  is  the 
pressure,  and  K  =      |_    grad  q  is  the  eddy  viscosity.    The  turbulent 
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1  I 

kinetic  energy,  E'  =    HL*  UA.     ,  is  assumed  to  have  the  form 

(2.4)  E'  =    C  |_2(grad  q)2 

where    C    is  an  unknown  constant.    It  is  assumed  that  the  rate  of 
dissipation  of  E'  is  proportional  to    E'. 

A  flux  of  turbulent  kinetic  energy  into  any  volume  element  is  assumed 
to  occur  when  the  quantity     U       varies  over  the  element.    Thus,  the 
rate  of  inflow  of    E1    is 


(2.5)  <=< 


(grad     L      ,     n)dS 


where    n    is  the  normal  to  any  surface  and    °^_      is  an  unknown  constant 
depending  on  the  turbulent  intensity.    The  turbulent  kinetic  energy  is 
constant  at  any  point,  and  the  above  assumptions  then  imply  that 

(2.6)  J   )    (grad    |_     ,  n)dS    +    B     1       E'  clV=0, 

where    B    is  an  unknown  constant.    Applying  the  divergence  theorem, 
considering    V    to  be  an  arbitrary  volume,  and  using  equation  (2.4), 
equation  (2.6)  becomes 

(2.7)  V2  L         +    B    L2(grad  q)2    =    0. 

Equations  (2.3)  and  (2.7),  together  with  the  continuity  equation, 
form  a  closed  set  for  the  mean  velocities  and  the  quantity    |_      .     The 
constant    B    must  be  determined  experimentally. 


For  flow  between  two  parallel  walls  as  shown  in  Figure  1,  the 
continuity  equation  is  trivial,  and  equations  (2.3)  and  (2.7)  become 


(2.8) 


-grad  P    + 


i.    |"L? 

dy  L 


{- 


V 


dy 


-) 


-    0 


and 


(2.9) 


dy 


+    B    LV   v     )2       =0 


dy 


Figure  1 


y=  1 


i 


y=  -1 


Defining    B    =  A    for  the  zero-order  flow  case  and  setting 
2 

equation  (2.9)  becomes 


=  |  grad  P  |  , 


(2.10) 


dy' 


J 4  A 


L 


(d 


T7 


dy 


0, 


and  equation  (2.9),   after  integrating,   becomes 


(2.11) 


hy 


2,H    V 


L> 


dy 


Equations  (2.10)  and  (2.11)  give 


(2.12) 


dy2 


-fhy. 
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which  has  the  boundary  conditions: 


d 


L 


dy 


L_  o     i£ 


=    0    and 


L 


L 


y  =  0  y  =  1 

where    i o    is  a  quantity  analogous  to  the  roughness  length  in  mixing- 
length  theory.    Integrating  equation  (2.12)  twice  and  applying  these 
boundary  conditions  results  in 


(2.13) 


L  -"ftV-a  +  lo. 


u 


Next,   solving  equation  (2.11)  for  d  \J       and  substituting  for 

dy 


L 


gives 


(2.14) 


TJ 


dy 


7h7 


L 


+  ^d  -y3) 

O  1^ 


By  defining 


=    l    +    12     L 

Ah 


and 


z    =    -*- 


VT 


equation  (2.14)  simplifies  to 


(2.15) 


12 


-1 


A/hF 


iTJ 


JT  dz 


1  -  z3 


Integrating  equation  (2.15),  applying  the  no-slip  condition  at  z 
and  returning  to  the    y  coordinate  gives 


=    -    N/~b7 
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(2.16) 


U(y)  = 


A/hF 


tanh"1 (I)     -    tanh*"1/  v3/2  , 


The  mean  velocity  at  mid-channel  is 


U(y) 


=      1 


y-  0 


which,  when  applied  to  equation  (2.16),  results  in  the  expression 


(2.17) 


aThf  tanh   m  '  x 


Equation  (2.16)  therefore  becomes 
(2.18) 


tanh^y^ 

U  (y)  =   i  -  VF 


ih_1(I 


L 


The  expression  (2.13)  for    L_      can  be  simplified  by  recalling  that 
von  Karman's  constant,     #V      ,  gives  the  rate  of  growth  of     I 


at  the 


wall.    Therefore, 


(2.19) 


where 


K  = 


L 


dy 


y-  -  1 


K« 


it  is  found  that 
(2.20) 
(2.21) 
(2.22) 


0.4  is  the  accepted  value.    Then,  from  equations  (2.13) 


Ah    =    1.6 

b      =    1    +    7.5  L 
.3 


L- 


b    -    v- 
7.5 
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Following  Benjamin  (1959),  the  equations  of  motion  in  terms  of    \ 
ind    ^7 


are 


(2.30a)  J<U,C^      "    t4  <^„„    )    +    l/2Tt(/^2    +(^2)  =  -P^     +  f, 

b) 


-^  ^  y,T  "^      -v^'-     -j     -^ 


J(-Kty   +  ^^  >  +  ^y^2  +Pt?2)  =  -i\  +  £ 


where  the  subscripts  denote  partial  differentiation. 

The  stream  function  is  composed  of  a  zero-order  term  plus  a 
periodic  perturbation.    When  there  is  no  wave,  the  stream  function  is 

'1r 


(2.31) 


H 


JJ(y[) 


d 


1 


where    "*|     +    y  =  1 . 

Vy      (V7)  is  the  velocity  profile  from  equation  (2.18)  expressed  in 
curvilinear  coordinates 


(2.32) 


TJ"i' 


tanh' 


1  - 


When  the  water  surface  is  disturbed  by  a  stationary  wave,  the  stream 
function  becomes 

(2.33)  ^     (fe,>})=U>        +aJF(>|)+  U  (7j)  ~c    je-^le11^ 

where  F(yi)  is  to  be  determined. 

The  mixing  length  is  expressed  in  a  form  similar  to  that  for  \y     : 

(2.34)  J[        =    [_(>f)     +    aHi(rl)eik^ 
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where     |_    (7f)  is  given  above  and  *X. ,  (  V»  )  is  to  be  determined. 
Similarly, 

(2.35)  B    =  \   +    a^T    oik) 

where    A    is  given  above  and    <j       is  an  unknown  parameter  to  be 
determined  empirically.    The  velocity  components  parallel  to    e       and   7? 
are  given  by 

(2.36a)       u    =Jl7Vn    =      U   -c  +  a  (f'CH)  +      U  (>j  )e~k>0-  elkJ? 

b)       v       *    2r         =    -ilea  J*  F(T,)  ♦    [TjT(IJ)  "  e]  e"^  }  eik*  . 


To  formulate  the  differential  equation  in  -JL,     to  be  solved, 
substitution  is  made  for  JC   ,     B    and  q  in  equation  (2.7)  expressed  in 
curvilinear  coordinates: 

(2.37)  vi        +    B(^_)2    =    0.. 

Simplification  by  elimination  of  all  terms  of  0(u2)  and  above  yields 

(2.38)  J[      +     JJ((AE      k2)    +    AHF"  +      ^D    + 

Ae~kyl    (G  -  kD)    =    0 


where 


*2    I     2 


d  -  U 2  L 
E  =  IP  L  . 


and 
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Similarly,  when  equations  (2.30)  are  expanded,  cross  differentiated 
to  eliminate  pressure  terms,  and  linearized  in    a,  one  obtains 


(2.39)  ik_ 

2 


i    (    \J     '  c)(F"  -  k2F)  -    "[J*  F     \    sgn(  ~[T  )    = 


(^—T    +    k2)(HF"    +    eJ?,)    + 
d7| 

e"kY(   (G"     -    2kG'    +    2k2G) 


where 


g=u'U"L2 

H =  rj'L2 

and  where    U       is  positive  for  the  lower  half  of  the  channel.     Note  that 
a  singularity  occurs  at  mid-channel  (    \J     =0).    Equations  (2.38)  and 

(2.39)  are  to  be  solved  for    F    and    j(try    . 

The  boundary  conditions  are  as  follows.    At  the  interface  where 
~*\     ~  0,        \j    (0)  =  0  to  satisfy  the  no-slip  condition,  and  v  =  0  because 
the  wave  is  stationary.     Therefore,  from  equation  (2.36b): 

(2.40)  F(0)    =    c. 

The  velocities  in  the  water  must  take  on  negative  values  in  order  that 
the  shearing  stress  across  the  interface  be  continuous.    The  variable 
component  of  the  tangential  water  velocity  associated  with  an  irrotational 
wave  is  of  the  form 


(2.41)  u'    =    a    x/T     /3H 


n 


17, 


2 


e 


ik 
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The  total  tangential  water  velocity  is 


(2.42)  u    -    -c    -    a/T  /sK_     eik£    . 

\J  O 


I 

o 


Defining  the  Froude  number 

2 

gH" 

tangential  air  velocity  is 


TT  2 

(2.43)  H~      =     % °     ,  and  invoking  the  no-slip  condition,  the 

n  h 


(2.44)  u    =    -c    -    a    l-^-     elk)        at     >7     =0,  or 


(2.45)  u    =    -c    -    a 


0  ^)     «  ??  =  o 


where  p    =  /~|p —  .    From  equation  (2.36a)  the  second  lower  boundary 
condition  then  becomes 


(2.46)  F'(0)    =    -     U'(0)    +     ft     . 

The  perturbation  in  Jc      at  the  water  surface,  JCt   (0) ,  must  also  be 
fixed.    This  must  be  regarded  (along  with  cT     )  as  an  unknown  parameter 
in  the  problem,  to  be  adjusted  empirically. 

The  perturbed  flow  must  die  away  at  some  large  distance  (relative 
to  k     )  from  the  wavy  surface.    Thus,  as  V?     becomes  large, 

(2.47)  Ft>|)    >-     (c    -U    )e"k7l 

(2.48)  F'(*W  ) H  -  \J*  e~k7i 

(2.49)  Jc,(>|) H  0. 
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2.3    NORMAL  PRESSURE,  SHEAR  STRESS  AND  WAVE  GROWTH 

The  pressure  on  the  wavy  surface  is  the  sum  of  that  due  to  normal 
atmospheric  turbulence  and  that  generated  by  air  moving  over  the 
undulating  surface.    The  latter  pressure  predominates  in  the  generation 
mechanism  considered  here  and  is  expressed  as 

(2.50)  P(0)    =    a(%     +    icf    )eikj? 

where  ^X/     and     o      are  to  be  determined.     Differentiating  equation 

(2.50)  with  respect  to    |     ,   substituting  the  result  into  equation  (2.30a), 
and  simplifying  yields 

(2.51)  P(0)    =       [THf    "    (    U     "    c)F'     -  ^     ^—     (HF" 

+    E^,)    +    (G'     "    Gk)    1 

J    -Y[  =  0 

which  gives  ^    and    o 

In  the  same  manner,  the  shear  stress  at  the  wall  is  expressed  as 

(2.52)  7T  (0)    =    a(6         +    i    £      )eik  Jf 

where    £        and     £,      are  to  be  determined.    To  evaluate  the  shear 
stress,  equation  (2.2)  is  expanded  in  terms  of    y      '• 

(2.53)  T    =    ^V  (^      +  ln  iyn  (y^). 

Substituting  for  JL    ,  ],      if/      ,  or  their  derivatives  where  required 
obtains 
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(2-54)  T  (0)    =    (    2    |_U"[|_U"         +     L  <F"  +   V"  >]) 


which  gives   £       and     £.     .    The  shear  stress  seems  to  be  relatively 
unimportant  to  wave  growth,  however,  and  will  not  be  considered  herein. 

To  compute  the  wave  growth  rate  the  normal  pressure  in  phase  with 
the  wave  slope  is  required.    Expansion  of  equation  (2.50)  and  comparison 
with  equation  (2.25)  for  the  wave  form  show  that  a  positive    o      results 
in  wave  growth.    The  average  energy  input  due  to  normal  pressure  is 

r* 

(2.55)  E    =    -   A       I         P^1-      d£ 

'o 

or,  in  dimensional  terms, 

~         a2  T72X 

(2.56)  E    =    |jjlO  fJJo     <*, 

2 
where     CO         =  gk  is  the  wave  frequency. 

The  surface  energy  of  the  wave  is 


(2.5  7)  E    =  1/2     O        ga2 


where     0         is  the  water  density ,   so  that 

(2-58)  I   -      J*,    ga^     . 

Equating  equations  (2.56)  and  (2.58)  and  solving  for  %$-  gives 


dT 


a  TT  2   r 

^•by>  dt  P  2gH 
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Integrating  produces  the  familiar  exponential  wave  growth  equation: 


(2.60) 


ao  exP    1       Y^  "  O  ^ 


^S 


-\ 


*"  J 

The  rate  of  energy  input  due  to  surface  shear  is 


r 


(2.61) 


1 


z% 


u     L      d\ 


/ 


or 


(2.62) 


E    =    a2 


£  -¥-^6. 


This  energy  is  expected  to  go,  at  least  initially,  into  a  surface  shear 
flow. 
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3.     SOLUTION  OF  THE  FIRST-ORDER  PROBLEM 

A  solution  for  the  function    F    and  its  derivatives  must  now  be 
found  so  that  values  for  normal  pressure,   shear  stress  and  wave  growth 
rate  may  be  determined.    There  are  several  numerical  integration  methods 
available  for  this  purpose,  having  the  common  property  of  starting  at 
one  boundary  and  ending  at  the  other.     However,  the  problem  as 
formulated  has  conditions  at  both  boundaries,  and  a  method  of  super- 
position of  solutions  must  be  used. 

Equations  (2.38)  and  (2.39)  can  be  rewritten  as  six  linear  first- 
order  ordinary  differential  equations.    That  is  letting: 


(3.1a) 

F     =Yl 

b) 

r  =y2 

c) 

F"   =  Y 

3 

d) 

F"    =  Y 

e) 

£=*, 

f) 

X  -  \ 

gives  six  ordinary  differential  equations: 


(3.2a) 

Yi   =  Y2 

b) 

Y*    =  Yo 
2          3 

c) 

K    =Y4 
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r 


d)    y;=hj 


2 


(U     -  c)  (Y3  -  k2Yx)  -     TJ"  Yx  -    2H'Y 

,2_      .    „£    _„„.„        „..     .,.2. 


-     (H"    +    k^H)Y0     -    E>L,      -  2E'YC  -  (E"    +  k*E)Yc 

O  DO 


-k 

-    e 


^     (2Gk2  -  2G'k    +    G")   f 


and 


e)  Y'    =Y_ 

5  6 

f)  Y'    =  -Y,-  (AE    -    k2)     -    AHYo    -        0    D 

6  o  6 

-Ae*"k?*    (G    -    kD). 

There  are  now  six  equations,  each  requiring  an  initial  input  at  the  lower 
channel  boundary.     Linearly  independent  solutions  are  found  by  starting 
repeatedly  at  the  lower  boundary  with  different  initial  values,  Y(0) .     The 
final  solution  is  obtained  by  superposing  the  independent  solutions  to 
fit  the  required  conditions  at  the  upper  boundary. 

The  independent  solutions,  denoted  by  the  superscript    j,  have 
initial  values  given  by 


(3.3)  Y^+l(o)    =       j     j     at    T(     =    0, 


where  is  the  Kronecker  delta.     The  solution  to  the  homogeneous 

form  of  equations  (3.2)  can  be  written  as  a  sum  of  independent  solutions 

to  the  homogeneous  equation: 

6 
<3-4)  Yi    =       y  h    Y- 

r-i       J     1 
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where  the    J's    are  complex  constants.    Since 

F(0)    =    F'(0)    =    Xx  (0)    -    0 
for  the  homogeneous  system,  the  constants  J-^,  J?,  and    J_    in  equation 
(3.4)  must  equal  zero.    The  complete  solution  is  achieved  by  adding  to 
the  homogeneous  solution  a  particular  solution  of  equation  (3.2).    For 
this  particular  solution,  the  initial  conditions  are 


(3.5a) 

Y^c 

b) 

Yp  -  $  -  -g 

c) 

YP    =    0 
3 

d) 
e) 

yp   =    0 
4 

Y^    =    i,(0) 

f) 

Yp    =    0. 

Then  the  complete  solution,   satisfying  all  the  boundary  conditions  at 

y\  =  o, 


IS 


(3.6)  Y.(>J)    =    T3^C|)    +    T4Y4(>()    +    T6Y6(^)    +     ^C?)' 

where  the  J's    are  to  be  determined. 

The  position  of  the  upper  boundary,  where  conditions  (2.47),    (2.48), 

and  (2.49)  hold,  is  not  known,  and  must  be  determined  by  a  convergence 

criterion.    The  solution  is  determined  by  applying  the  upper  boundary 

condition  at  various  increasing  values  of  "^     =    d ,   solving  the  resulting 

equations  (3.6)  evaluated  at    d    for  the    J's,     and  evaluating  the  wall 

pressure  each  time  using  equation  (2.51).    When  the  wall  pressure 
determined  in  this  manner  converges,  the  procedure  is  stopped. 
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4.     THE  COMPUTER  SOLUTION 

The  solution  for  the  function    F    and  its  derivatives  has  been 
obtained  by  solving  the  problem  as  outlined  in  section  3  on  an  IBM  360- 
67  computer  using  complex  double-precision  arithmetic.     Major  sub- 
routines are  incorporated  in  the  program  to  calculate  the  velocity  profile, 
perform  the  integration  using  the  Runge-Kutta  Gill  fourth-order  method, 
compute  the  shear  stress,  calculate  the  wave  growth  rate  and  evaluate 
the  stream  function.    Appendix  I  contains  the  format  of  the  program  and 
a  set  of  results  obtained. 

The  accuracy  of  the  Runge-Kutta  Gill  integration  method  depends 
on  the  size  of  the  step  used  in  the  iteration.    The  value  10       was 
selected  in  order  to  determine  whether  the  surface  pressure  converged 
on  a  particular  value.     Because  surface  pressure  did  not  converge  for 
this  step  size,  it  was  decreased  until  convergence  was  achieved.    A 
step  size  of  10       was  necessary  to  obtain  convergence  of  the  wall 
pressure,  and  the  program  running  time  was  about  112  minutes.    A  final 
optimum  step  size  of  10   °  for    >]      less  than  0.3  and  2.5  X  10~5  for    V 
greater  than  0.3  was  then  adopted.    This  step  size  gives  convergence 
on  a  result  essentially  the  same  as  that  for  10       and  a  running  time 
of  65  minutes . 
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RESULTS  AND  CONCLUSIONS 


The  zero-order  velocity  profile  is  given  by  equation  (2.18)  which 

has    I (0)  as  a  parameter.    A  plot  of     U    (y)  for  several  values  of   I — (0) 

appears  in  Figure  3  .     From  the  figure  it  can  be  seen  that  as  the  parameter 

| (0)  decreases,  the  shape  of  the  profile  changes  from  nearly  linear  to 

almost  a  step  function.    A  representative  logarithmic  profile  is  also 
shown  in  the  same  figure  so  that  the  two  profile  shapes  may  be  easily 

compared. 

1.0     f 


0.8     - 


0.6 


0.4     . 


0.2     - 


1.0 


Figure  3 . 
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The  remainder  of  the  analysis  is  concerned  with  the  determination 
of  suitable  values  for  the  parameters     [__  (0) ,     Jl,   (0) ,  and     o        .     For 
the  channel  under  consideration: 
H    =    10  meters 
\J0=    10  meters/second 
which  results  in  a  Froude  number  of  1.018.     For  the  wave  number, 

k    =    15, 
the  instability  mechanism  considered  is  important.    These  values  for 
the  Froude  number  and  the  wave  number  have  been  used  throughout  the 
analysis. 

The  following  parameter  values  were  first  used: 
|__(0)    =    2"22  X  10~2 
J2,(0)     =    0 
**         =    10. 
The  wall  pressure  failed  to  converge  for  these  values.     Convergence 
was  assumed  to  have  resulted  if  wall  pressure  values  were  nearly 
constant  over  a  significant  range  of    d  (say  0.1)  or  reached  a  minimum 
value  at  some  point  in  the  channel.     (This  latter  criterion  is  not 
completely  defensible,  but  the  singularity  in  equation  (3.  2d)  at  77     =  1, 
itself  a  peculiarity  of  the  Zagustin  assumptions,  is  probably  responsible 
for  the  erratic  behavior  of  the  wall  pressure  as  77— KL.)       0      was  then 
varied  from  -1000  through  zero  to  +1000  and  given  imaginary  values  as 
well.    The  wall  pressure  changed  only  in  the  fourth  or  fifth  decimal 
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which  was  considered  insignificant,  although   ^-i 

Thus,  the  solution  seems  insensitive  to     0 

New  values  for  the  parameters  were  then  used: 

,-5 


changed  markedly, 


L  (0)    =    5X1 
2t  (0)    =    0 
&         =    10. 


The  wall  pressure  converged  for  these  values,  but  the  value  was  too 
high  in  magnitude  and  of  the  wrong  sign.    A  value  nearer  to  5  (Miles, 
195  7)  was  considered  appropriate. 


With 


L  (o) 

=    5  X    10 

&  (o) 

=    1 

y 

=    10 

-5 


the  resulting  wall  pressure  was  positive  but  still  too  high  in  magnitude. 
Figure  4  summarizes  the  results  of  varying       | (0)  and    ->L,  (0). 


1    d  (position  where  upper  boundary 
1 .  0  conditions  are  applied) 


-«* 


l_(0)  =  2.22  X  10 
J?,(0)  =  0 

y   =iq 

1  win  "' '  > 

I (0)  =  5  X  10 

J2i(0)  =  o 

=  10 


•— < — 

-io' 


-10J         -10^         -io 

Wall  Pressure, 
Figure  4 . 
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It  appeared  that  the  wall  pressure  was  more  directly  dependent  on 
Jit  (0)  than  on  any  other  parameter.    Subsequent  investigation  showed 
that  for  any  >?      the  wall  pressure  was  a  linear  function  of  the  real  or 
imaginary  part  of   *X.X  (0)  providing  that  the  other  component  was  held 


constant.     For 


Leo 

2,(0) 


5  X  10"5 
0.65788 
10 


the  resulting  minimum  wall  pressure  was  3.83  at    d  =  0.48.     Figure  5 
shows  the  resulting  wall  pressure  in  this  case  versus    d. 


1.0 

0.8 

0.6 

0.4 

0.2 


10 


5  X  10~5 
0.65788  +0i 


=  10 


-w* 


J h 


102     103      104 

Wall  Pressure,  O 
Figure  5 . 


0      was  then  given  a  value: 

S*    =  1000. 
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-^.(0) 

Out-of-phase 

In-phase 

Wall  Pressure 

,i 

Shear  Stres 

0.0     4- 

Oi 

-6495.3 

-29.78 

0.2     + 

Oi 

-4519.4 

-14.89 

0.4      + 

Oi 

-2543.6 

0.001 

No  significant  change  occurred  in  the  wall  pressure.     Convergence,  as 

defined  above,  always  occurred  near    d  =  0.5. 

Table  1  summarizes  the  values  of  the  wall  pressure,   shear  stress, 

and  wave  growth  rate  versus  Jl,  (0)  at    d  =  0.5  for  parameter  values: 

L(0)    =    5  X  10"5 

r       =    10. 

Wave  Growth 
ar  Stress, 6.  Rate ,  j^  iE_  X ^ 

-1.61 

-1.12 

-0.63 

0.6     +     Oi  -567.7  14.87  -0.141 

0.349 

0.838 

0.001 

-1.61 

-1.58 

-1.55 

-1.52 

-1.49 

0.013 


Table  1.         Wall  pressure,  shear  stress,  and  wave  growth  rate 
as  a  function  of  the  perturbation  surface  roughness 


0.8     + 

Oi 

1408.2 

29.76 

1.0     + 

Oi 

3384.0 

44.64 

0.65788+0i 

3.83 

19.19 

0.0     4- 

0.21 

-6495.3 

-29.78 

0.0     + 

0.41 

-6401.8 

-29.78 

0.0      + 

0.61 

-6308.3 

-29.78 

0.0     + 

0.81 

-6214.8 

-29.78 

0.0     +- 

1.01 

-6221.3 

-29.78 

0.65788+0. li 

50.9 

-29.78 
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It  can  be  seen  that  imaginary  values  for   ^-x  (0)  have  no  effect  on  the 
shear  stress.    The  required  value  for    J.,   (0)  is  0.65  788  +  Oi  for  a  wall 
pressure  of  3.83,  but  the  shear  stress  seems  high.     (However,  the  shear 
stress  value  is  a  first-order  contribution  to  the  total  shear  stress  and  is 
reduced  by  the  fractional  value  of  the  wave  amplitude.)    Exponential 
wave  growth  rates,  due  to  the  wall  pressure  alone,  appear  reasonable. 
To  reduce  the  shear  stress  to  a  more  acceptable  value  in  relation  to  the 
normal  pressure,  parameter  values 


L  (0)    =    5X10   5 


X,  (0)    =    0.4  +  5.45i 

8*      =  10 

were  used,  and  give  at    77     =    0.5 

Normal  pressure  =    5.0 

Shear  stress  =    0.001 

Growth  rate  =    0.001. 

However,  the  wall  pressure  failed  to  converge,  and  had  the  form  shown 
in  Figure  6. 

Equation  (2.33)  for  the  stream  function  has  been  evaluated  with 
the  values  of    F    calculated  from  the  constants    J    determined  at  d  =  0.48 
where  the  wall  pressure  was  found  to  converge.     Parameter  values 


were 


|_(0)    =    5X10   5 
Jl,(0)    =    0.65788 
V         =    10. 
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Figure  7  shows  a  plot  of  the  streamlines  with  the  'cats  eye1  present 
and  centered  about  30  degrees  ahead  of  the  wave  crest.    The  'cats  eye' 
or  closed  loop  is  indicative  of  the  wave  induced  perturbations  in  the 
flow  and  characterizes  the  region  of  turbulence  (Phillips,    1966).     Higher 
in  the  channel  the  streamlines  become  less  wavy  and  approach  straight 
lines.    Table  2  shows  values  of  the  stream  function  near  the  critical 

* 

layer. 

In  view  of  the  foregoing  limited  results,  the  following  tentative 
conclusions  can  be  made: 

1)  The  inverse  hyperbolic  tangent  velocity  profile,  a  result  of 

the  theory  of  the  model,  can  be  fitted  to  almost  any  known  profile  type. 
The  profile  which  has  a  large  value  for    — VJ  near  the  critical  level 

and  for  which  the  parameter    {_  (0)  has  the  value    5X10         produces 
surface  pressures  which  can  be  made  to  converge  on  the  desired  value 
by  varying  the  wall  perturbation  in  ~M.  . 

2)  Wall  pressure  is  extremely  sensitive  to  very  small  changes  in 
-^i(O)  and  insensitive  to  moderate  changes  in    o       . 

3)Reasonable  wall  pressure  (and  wave  growth  rates)  can  be  obtained 
with  an  appropriate  choice  of   JCX  (0).    This  indicates  that  small,  unequal 
variations  in  the  surface  roughness,   such  as  would  be  caused  by 
parasitic  capillaries,  are  very  important  in  the  wave-generation  process. 

4)    Shear  is  independent  of  imaginary  values  of    X,  (0).    Therefore, 

to  vary  the  shear  stress,    0\jfc-    (   ^  ,  (0)),     0      or    I (0)  must  be  varied 

individually  or  in  combination  as  indicated  by  further  analyis  of  the  model, 
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Based  on  the  above  results,  the  following  suggestions  are  made: 

1)  I (0)  should  be  increased  slightly  to 

L(o)  =  io"4 

in  order  to  reduce  the  shear  stress  relative  to  the  wall  pressure.    One 
computer  run  indicates  that  the  shear  stress  is  halved  for  this  value. 

2)  Additional  wave  numbers  should  be  investigated,   say,  in  the 
range  5   C     k    C    20,  and  their  effect  on  wall  pressure  and  shear  stress 
noted. 

3)  The  model  should  be  checked  using  data  from  Stanton  and 
Motzfeld  as  quoted  in  Ursel?  (1956),  and  with  recently  obtained  laboratory- 
data  (see,  for  example,  Zagustin,  et  al.  ,   1966). 

In  summary,  the  model  appears  to  be  fairly  well  suited  to  a 
theoretical  study  of  air-stream  turbulence  and  wave  generation. 
Unfortunately,  the  long  time  required  for  a  computer  run  prevented  the 
calculation  of  further  data  using  appropriate  parametric  values. 
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k 

ChO 0_A (h_2 CLJ 0.4 

0.0  0.0  0.0  0.0  0.0  0.0 

0.002  -1.522  -1.256  0.508  0.434  0.216 

0.005  -0.683  -0.571  -0.227  0.218  0.594 

0.008  -0.299  -0.253  -0.080  0.154  0.359 

0.010  -0.138  -0.116  -0.007  0.145  0.285 

0.015  0.110  0.101  0.129  0.183  0.243 

0.020  0.244  0.229  0.238  0.269  0.309 

0.025  0.325  0.315  0.336  0.381  0.432 

0.035  0.424  0.440  0.525  0.646  0.758 

0.045  0.506  0.556  0.722  0.940  1.128 


H 

0.5 

0.6 

0.7 

0.8 

0.9 

1 

0.0 

0.0 

0.0 

0.0 

0.0 

0.0 

0.002 

1.533 

1.266 

0.518 

0.426 

-1.522 

0.005 

0.757 

0.645 

0.300 

-0.145 

-0.520 

0.008 

0.456 

0.410 

0.236 

0.003 

-0.202 

0.010 

0.357 

0.336 

0.228 

0.075 

-0.065 

0.015 

0.284 

0.293 

0.265 

0.211 

0.152 

0.020 

0.343 

0.359 

0.349 

0.319 

0.278 

0.025 

0.470 

0.480 

0.459 

0.414 

0.363 

0.035 

0.817 

0.801 

0.716 

0.594 

0.482 

0.045 

1.212 

1.162 

0.996 

0.778 

0.591 

Table  2.       Values  of  the  stream  function  near  the  critical  layer 
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